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Dendritic morphology constrains brain activity, as it determines first which neuronal 
circuits are possible and second which dendritic computations can be performed over 
a neuron's inputs. It is known that a range of chemical cues can influence the final shape 
of dendrites during development. Here, we investigate the extent to which self-referential 
influences, cues generated by the neuron itself, might influence morphology. To this end, 
we developed a phenomenological model and algorithm to generate virtual morphologies, 
which are then compared to experimentally reconstructed morphologies. In the model, 
branching probability follows a Galton-Watson process, while the geometry is determined 
by "homotypic forces" exerting influence on the direction of random growth in a 
constrained space. We model three such homotypic forces, namely an inertial force 
based on membrane stiffness, a soma-oriented tropism, and a force of self-avoidance, as 
directional biases in the growth algorithm. With computer simulations we explored how 
each bias shapes neuronal morphologies. We show that based on these principles, we can 
generate realistic morphologies of several distinct neuronal types. We discuss the extent 
to which homotypic forces might influence real dendritic morphologies, and speculate 
about the influence of other environmental cues on neuronal shape and circuitry. 
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1. INTRODUCTION 

Dendrites are beautiful arbors sprouting from the cell bodies 
of neurons. The shape of dendrites is of great importance to 
the nervous system for two interrelated reasons. First, a neu- 
ron's dendrites receive inputs from other neurons. Since dendritic 
morphologies define which spatial domain can be reached, they 
govern the connectivity of the local circuit. Second, neuronal 
computation takes place along the entire dendritic path between 
the site of the input and the site of action potential genera- 
tion in the axon (Koch and Segev, 2000; London and Hausser, 
2005; Torben-Nielsen and Stiefel, 2010). Dendrites exist in a wide 
range of sizes and shapes and their morphologies have historically 
been used to classify neuronal types (Hillman, 1979; Migliore and 
Shepherd, 2005). Moreover, dendritic trees of one neuron type 
can display a large variation (Soltesz, 2005), and in some cases 
variations in the morphology of the dendritic tree is indicative of 
particular neurological disorders (Irwin et al., 2000; Kaufmann 
and Moser, 2000). Consequently, the morphology of dendrites 
has been, and remains, at the core of many studies. 

To assign neuronal morphologies and their traits to distinct 
neuronal types, invariant descriptors are required. To this end, 
distributions of properties are measured directly from digitally 
reconstructed neurons (for instance, the distribution of segment 
lengths). Parameterizations of these distributions (Ascoli and 
Krichmar, 2000; Samsonovich and Ascoli, 2005; Koene et al., 
2009) or the distributions themselves (Lindsay et al, 2007; 



Torben-Nielsen et al., 2008) then provide these invariant descrip- 
tors of morphology. However, despite considerable success, these 
approaches cannot account for all experimentally observed vari- 
ations in the data. Indeed, it is often impossible to construct 
an exhaustive statistical model because the experimental data 
does not contain sufficient samples, or is simply not coherent 
(Torben-Nielsen et al, 2008). 

An interpretation thereof is that a single dendritic morphol- 
ogy is the result of interactions with the environment in which 
it grew. These interactions are de facto not in the reconstructed 
data, which is a snapshot of one of their possible outcomes. 
Experimental evidence backs this interpretation as it has been 
shown that dendrites are shaped by chemical cues in the environ- 
ment (Scott and Luo, 2001; Grueber and Sagasti, 2010; Jan and 
Jan, 2010), and a specific class of these environmental interactions 
are termed homotypic or self-referential, that is, originating from 
the neuron itself. 

At least three self- referential forces derived from these cues 
are experimentally and theoretically described, namely self- 
avoidance (Grueber et al, 2005; Marks and Burke, 2007b; 
Grueber and Sagasti, 2010; Jan and Jan, 2010), soma-oriented 
tropism (Samsonovich and Ascoli, 2003; Marks and Burke, 
2007a), and mechanical stiffness, which provides a bias to the 
growth cone to favor straight motion (Condron and Zinn, 1997; 
Koene et al., 2009). In this work, we propose a novel method, 
extending the work of Samsonovich and Ascoli (2003) and 
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Marks and Burke (2007b), to model these self- referential inter- 
actions that shape neuronal morphologies. The method algo- 
rithmically generates dendrograms based on a Galton-Watson 
process and simultaneously adds geometry that is locally deter- 
mined solely by the three aforementioned directional biases. 
We show that based on these simple principles, virtual mor- 
phologies with high resemblance to neurons can be generated. 
Hence, we demonstrate in silico that self- referential cues are 
sufficient to shape dendritic morphologies realistically given 
an otherwise random growth process and are thereby capable 
of generating multiple isometric variants of a single recon- 
structed instance. Because self- referential cues can account 
for these properties, we propose their descriptions might be 
useful as an integral part of overall neuronal morphological 
descriptions. 

2. METHODS 

2.1. M0RPH0GENETIC ALGORITHM 

We present an integrated morphogenetic algorithm (Figure 1), 
in which a neuron's dendrogram and geometry are generated 
simultaneously, such that the dendrogram derives from a mod- 
ified Galton-Watson process, while the local geometry is dic- 
tated solely by the sum of self- referential growth directional 
biases. Neurons and their branches are further subject to cer- 
tain termination conditions. In addition, we include a random 



component of the morphogenetic process, a 3-dimensional (3-D) 
Gaussian distribution, from which all directions of growth are 
sampled. 

A simulation begins with a specific configuration of model 
parameters (Table 1; we include an exemplary configuration file 
for each simulation in the Supplementary Materials). The main 
parameters of the algorithm describe the strength, spatial gradi- 
ent, and extent of the local growth biases, and the branching and 
extension processes. A second subset of parameters describes the 
termination conditions for both growth of individual branches 
and growth of the neuron as a whole. A final subset includes aux- 
iliary parameters for initial conditions (e.g., soma surface area, 
number of stems, etc.). While termination and initial conditions 
are typically derived directly from experimental observations, 
parameters governing growth biases were chosen by hand to pro- 
duce morphologies that match other secondary measures such 
as space coverage and fractal dimension, as well as qualitative 
observations from experimental reconstructions. While both ter- 
mination and initial conditions had a potent effect on global 
properties of the generated cells' morphologies (size, total fiber 
length, etc.), they were less able than local homotypic growth 
biases to alter those cellular morphological traits that often define 
distinct types. These traits are further demonstrated by artificial 
morphologies generated when holding termination and initial 
condition constant (Figures 2, 3). 



Create soma at origin 

i 



M 



Read parameter file 



<«— | START ^ B 



Create n neurite stems with 
angles <minInitialStemAngle 



Extension Of Non-terminated Dendrite Ends 



Application of Forces 




Apply forces to new front(s): 

• Inertial 

• Soma-oriented tropism 

• Self-avoidance 



-[ Grow Segment(s) ] 



Compute Forward Direction 



Center 3-D Gaussian on front origin 



Offset to apply inertial force 



Compute Soma Direction 



Offset to apply soma-oriented tropism I 



Compute Vector Sum of Neighbors 



T 



Offset to apply self-avoidance 



Sample New Direction 



Output neuron swc file ]- 



FIGURE 1 | Schematic of the morphogenetic algorithm. The 

phenomenological algorithm relies on a Galton-Watson process to 
create a topology while the geometry results from applied forces from 



the environment. (A) Main algorithm to generate dendritic 
morphologies. (B) Procedures to sample angles biased by self-referential 
forces. 
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Table 1 | Parameters of the morphogenetic algorithm and their explanation. 





self_avoidance_force 


Maximum strength of self-avoidance bias: 0 indicates no self-avoidance, positive/negative (unbounded) values 




introduce a bias 


self avoidance decay 


Self-avoidance force decays as a function of distance to this power 


soma_tropic_f orce 


Maximum strength of soma-tropic bias: 0 indicates no soma-tropic bias, positive/negative (unbounded) values 




introduce a bias 


s oma_t r op i c_de c ay 


Soma-tropic force decays as a function of distance to this power 


inert ial_f orce 


Maximum strength of forward growth bias due to stiffness: 0 indicates no bias, positive (unbounded) values introduce 
a bias 


br anch_pr obab i 1 i t y 


Global probability of bifurcation over one unit length of front. (Probability of branching at a front is equal to the 




probability of branching at least once over total front length; see section 2.3) 


f r ont _ext ens i on 


Parameter for varying the length of a front's segments: 0 indicates "fixed" and otherwise "surface area dependent." 




With positive surface area dependence, front segments becomes longer as their radius decreases; with negative, 




shorter 




bounding_shape 


Parameterized surfaces at which growing branches terminate 


intersect ion_proximity 


Minimum distance between an active front and any other segment 


max_f iber_length 


Maximum total fiber length in the neuron, at which all branches terminate 


max_bi f urcat ions 


Maximum number of bifurcations in the neuron, at which all branches terminate 


min_radius 


Minimum segment radius, at which a branch terminates 


AUXILIARY PARAMETERS 




nbr_s terns 


Number of soma branches 


start_radius 


Starting radius of each dendrite 


taper_rate 


Decay constant of segment radius over sequential fronts 


ralls_ratio 


Ratio of the sum of initial branch segment radii to their terminal parent branch segment radius 


angle_limit 


Lower an upper limit on the bifurcation and stem angles. Limits on bifurcation angle is rarely used but insures biological 




plausibility. Minimum of the uniform distribution of angles between inertial force vectors resulting from a bifurcation 


flatness 


Scaling factor for the z-dimension of 3-D Gaussian distributions from which growth directions are sampled: 0 




generates a 2-D neuron, 1 generates a neuron with unconstrained 3-D morphology. (For flatness <0.5, bifurcation 




angles are sampled in the xy plane.) 



2.2. ALGORITHM STEPS 

After parsing the parameter text file, the C++ program that 
implements our algorithm first instantiates the soma with 
the specified surface area. The required neurite stems (soma 
branches) are created by iteratively sampling a random direc- 
tion from the soma, calculating the angle to all other stems, 
and adaptively decreasing the minimum angle until an appro- 
priate set of angles, and number of stems, is achieved. The algo- 
rithm then proceeds by repeatedly considering non-terminated 
branches ("fronts") and extending each by either continuing or 
branching. Newly created branches are pushed onto the list of 
non-terminated branches (fronts). When continuing, the branch 
grows in a direction determined by adding a vector representing 
each directional bias to a sample from the 3-D Gaussian distribu- 
tion, thus modeling the sum of the three homotypic forces. 

When branching, two daughter branches are subjected to 
these biases, calculated separately for each. The algorithm revisits 
all non-terminated branches repeatedly until a neuron termina- 
tion condition is reached. Termination conditions for branches 
include, in order of importance: (1) minimum segment diameter 
reached, (2) close proximity to another branch (measured sur- 
face to surface), and (3) collision with a fixed neuron boundary. 



Termination conditions for all neuron growth included, in order 
of importance: (1) all neuron branches terminated individu- 
ally, (2) maximum total fiber length reached, and (3) maximum 
number of bifurcations reached. Finally, the morphology data is 
ordered and outputted according to the "swc" format (Cannon 
et al, 1998). 

2.3. IMPLEMENTATION DETAILS 

2.3. 1. The random process and force modeling 

The influence of each of the three homotypic forces is modeled as 
an offset to the mean of a 3-D Gaussian distribution from which 
a direction of growth is sampled (Figure 1). For each dimension, 
the algorithm draws a sample from a Gaussian distribution mul- 
tiplied by a constant standard deviation parameter for the neuron 
and summed with a variable offset for that dimension determined 
by the sum of local growth biases, where each unbiased sample, s, 
is calculated using Marsaglias polar method (Marsaglia and Bray, 
1964) as follows: 

/ -2ln(vf + v 2 2 ) 

s = ovi , (1) 

V ^ + y2 
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where a is the standard deviation and {vi, v{\ e U(— 1, 1)|0 < 
v i + v 2 < wnere ^( a > is me uniform distribution over the 
interval [a,b]. 

This process introduces stochasticity into our generative 
model of neuron morphology (Figure 2A). For certain mor- 
phologies (e.g., alpha motor neurons vs. Purkinje cells), the 
randomness of growth direction was increased by increasing the 
standard deviation of the neurons 3-D Gaussian distribution. 
Furthermore, because of this stochasticity, a single parameter set 
could be used to generate a very large population of unique mor- 
phologies by varying each instance's random seed. (The resulting 
neurons have similar morphological traits and measures but are 
not identical.) 

For each new segment added by continuing or branching at 
a front, a Gaussian distribution is centered on the front origin 
(see Equation 1). Next, a vector for the inertial force is calculated 
as a continuation of the previous segment, or as the direction of 
the randomly generated bifurcation angle. The magnitude of this 
vector is scaled by a user-specified bias-specific parameter, and 
the mean of the Gaussian displaced accordingly. A vector for the 
soma-tropic force is calculated similarly, representing the direc- 
tion from the soma to the front origin, and additionally scaled 
by a user-specified function of this distance. Finally, a vector 
sum of the directions from all existing neuron segments to the 
front origin, with each direction vector scaled by a second user- 
specified function of distance, is calculated. With this vector sum, 
the mean of the Gaussian is displaced a third time. The center of 
the Gaussian may then be far from its initial mean. 



Since only the direction is sampled from this Gaussian, the 
greater the displacement of the mean, the more likely the actual 
direction vector will be in the direction of the sum of the three 
forces, though this is never guaranteed since the Gaussian is of 
infinite extent. Finally, a segment is added in the sampled direc- 
tion, and its length assigned according to the user- specified rate 
of growth. By adjusting these three forces, we show that the 
stochastic process may be biased to generate a variety of neuron 
morphologies (Figure 2B). 

2.3.2. Rate of growth and effects on morphology 

Both continuing and branching at the end of a branch are condi- 
tioned by the rate of growth, which may be specified as a constant 
rate (in units of segment length per front, resulting in a pure 
Galton-Watson process), a rate which depends on the variable 
(tapered) radius (in units of segment surface area per front), or 
some intermediate rate. For variable growth rates, the Galton- 
Watson process is modified, such that branching occurs after a 
segment of growth according to the branching probability, con- 
ditioned by the segment's length. In this way, a single bifurcation 
point occurs at the end of each segment if the probability of at 
least one bifurcation point over that length exceeds a randomly 
generated number on the unit interval. Experimentally, it has 
been shown that in particular neuron types, the branching proba- 
bility is non-uniform across the dendritic tree (Burke et al., 1992; 
Nowakowski et al., 1992). We found that certain non-uniform 
branching probabilities can be better approximated by our mod- 
ified Galton-Watson process when we specify a non-uniform 
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FIGURE 2 | Homotypic forces can shape dendritic morphologies. All 

illustrations are 2-D projections of 3-D structures. (A) Branched structure 
resulting from a Galton-Watson branching process without homotypic 
forces resembles a random diffusion process. (B) Dendritic-zlike structures 



emerge when different homotypic growth biases are added to define the 
geometry. The influence of different levels of inertial, soma-tropic, and 
self-avoidance are shown. Structures are bounded by a cube with side 
dimensions of 50|xm. 
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FIGURE 3 | Dendritic morphological traits associated homotypic forces. 

2-D dendritic morphologies generated with different settings for the 
homotypic growth biases. The modeled strength of one force is varied while 
the other two forces are fixed to a low strength. Top row: Dendrites with 
distinct levels of inertial forces. Middle row: Dendrites with distinct levels of 
soma-tropic forces. Bottom row: Dendrites with distinct levels of 



self-avoidance forces. Morphologies generated with strong inertial forces 
show sparse dendrites, while dense arbors require strong soma-tropic or 
self-avoidance forces. Inertial forces and self-avoidance grow larger 
structures in the same surface (space). Self-avoidance covers the surface 
(space) most densely. Structures are bounded by a rectangle with 
dimensions 150 x 100 |xm. 



growth rate. In addition, because segments are straight, we found 
that certain neurons whose tortuosity decreases or increases with 
distance from the soma can also be approximated by specifying 
a non-uniform growth rate under constant influence from our 
three forces applied to segments of increasing length. 

2.3.3. 3-D vs. 2-D growth 

Certain neuron types, such as retinal ganglion cells, starburst 
amacrine cells, and Purkinje cells, exhibit morphologies which 
are flat, covering a surface rather than filling a space. Our sam- 
pling procedure was initially developed for 3-D growth, but 
because growth depends entirely on sampling a 3-D Gaussian 
distribution, we were able to modify the procedure in a straight- 
forward way in order to restrict growth to two dimensions by 
modifying the Gaussian distribution. More specifically, the sam- 
pling procedure was adjusted to include a user- specified scaling 
parameter applied to the standard deviation of the distribu- 
tions z- dimension. In this way, the direction of growth could be 
constrained in this dimension, or eliminated altogether (with a 
scaling factor of zero). In addition, we modified the procedure for 
determining the direction of inertial force vectors after branch- 
ing, such that for scaling factors <0.5, bifurcation angles were 
restricted to the xy plane. 

2.3.4. Morphometric analysis 

We analyzed the resulting morphologies' morphometric prop- 
erties. In keeping with the aims of this work, we were most 



concerned with type-specific traits, such as the dimension of 
a morphology, the total length, the average path length to the 
terminal tips, and, space-filling proxies such as the fractal dimen- 
sion and the contraction of branches. The fractal dimension is 
computed on a branch-by-branch basis by the Hausdorff fractal 
dimension (as described in Marks and Burke, 2007a): ((E)) D = t, 
with (E) being the average increase in Euclidean distance per 
segment length and D the actual fractal dimension. A straight 
line always has the same Euclidean excursion from the soma per 
unit of path length and will result in D = 1, while a completely 
random walk will cover the surface (in 2-D) and have D = 2. 
Contraction is a similar metric indicating the ratio between the 
Euclidean distance and path length to a given terminal tip. The 
total number of branch points and terminal tips are also analyzed. 
For direct visual comparison of the topology we use dendrograms, 
and measured the asymmetry, A, of the tree of degree n, by sum- 
ming over all previous branch points, i , each of which creates two 
subtrees with terminal segment counts r and s, according to van 
Pelt et al. (2001), 

1 n ~ l 

A= -VApfesi) (2) 

n — 1 ^ 

i=l 

where A P (r, s) = \r - s\/ (r + s - 2) for {r, s} # {1, 1}, and 0 
otherwise. Finally, a Sholl-like analysis was performed to inves- 
tigate the distribution (conditioned on branch order or path 
length) of selected morphometric properties (Sholl, 1953). 
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2.3.5. Parallel implementation and availability 

The compiled algorithm is very fast, with each neuron morphol- 
ogy presented here generated in <2s using a 2.33 GHz micro- 
processor. For neural tissue simulation, in which thousands to 
millions of neurons may be generated in a single defined tissue 
volume (Kozloski, 2011), we anticipated that serial neuron mor- 
phogenesis will be too costly, and have therefore mapped our 
algorithm onto the parallel architecture of IBM's Blue Gene/P, 
such that each core of this massively parallel machine indepen- 
dently generates a unique set of neuron morphologies. By varying 
the random seed used to initialize this process on each machine 
core, we can produce a large space-filling tissue represented by 
a large number of swc files written in parallel to the machine's 
file system. These can then be read into a neural tissue simu- 
lator (Kozloski and Wagner, 2011) for subsequent physiological 
modeling. The final software is experimental. IBM would like 
to create an active user community for its neural tissue sim- 
ulation tools. Readers are therefore encouraged to contact the 
corresponding author if interested in using the tool or in the 
source code. 

3. RESULTS 

3.1. MORPHOLOGICAL TRAITS CAUSED BY H0M0TYPIC FORCES 

The influence of self- referential forces are manifested at both 
bifurcation points and continuing dendritic sections. As a refer- 
ence, we generated a morphology in which no homotypic forces 
were present (Figure 2A). The resulting morphology resembles 
no neuron and instead represents a diffusion process (Brownian 
motion) modified by particle division and elimination according 
to our rules for branching and termination. We then added the 
three different self- referential forces at various levels of strength. 
The resulting morphologies are shown in Figure 2B. We observe 
that the introduction of homotypic forces unclutters the ran- 
dom morphology and brings structure to the shape. Moreover, 
the morphologies subjected to the homotypic forces resemble 
dendritic morphologies. 

Next, we systematically investigated the morphological traits 
associated with a particular force. For the sake of clarity, we 
present neurons generated solely with 2-D growth, since the 
resulting dendritic trees are easier to inspect visually. Figure 3 
summarizes these results. The strength parameter of each applied 
force is varied from 1 to 50 (arbitrary units; Table 1). For each 
force combination we generated 25 separate morphologies from 
different seeds to gather population statistics (although only one 
is depicted; Figure 3). 

At the lowest level studied for each bias, we applied all three at 
strength 1 to impose structure (see Figure 2A). Therefore, mor- 
phologies in the first column of Figure 3 are similar as they all use 
the same force parameters (but different seeds). As a population 
these "baseline" structures have a total length of 2970 /xm ± 800 
and 230 =b 60 bifurcation points. The fractal dimension based 
on the Hausdorff metric is not always defined (Cannon et al., 
1999), but the closely related contraction ratio was measured at 
0.79 ±0.05. 

Next, we investigated the effect of increasing one of the three 
self- referential forces. First, in setting a higher inertial force 
parameter of 5 (top row in Figure 3), the total length initially 



increases to 4352 |xm zb 800 and 327 =b 63 branch points. At 
the highest parameter setting of 50, the size decreased again 
to slightly above the base line measurement (i.e., measurement 
when all force parameters set to 1). Then, we investigated the 
influence of the soma-tropic force and an entirely different mor- 
phology appears. With the highest soma-tropic force parameter 
used, the total length increases to 6388 |xm zb 2299 and the mor- 
phology consists of 488 ±176 branch points. The contraction 
increases to 1, which means that the dendrites follow a straight 
path directly away from the soma. Given the drastic increase in 
length when modeling the soma-tropic force rather than the iner- 
tial force, we can say that neurons subjected to this force form 
denser dendrites. Finally, when investigating the influence of the 
self- avoidance force, we see that the dendrites only become even 
denser. More specifically, with the highest tested self- avoidance 
parameter, the morphologies have a total length of 1 1, 293 |xm ± 
1181 and 808 ±144 branch points. The contraction at this level 
is 0.94 ±0.2. 

We thus observe that strong soma-tropism and strong self- 
avoidance results in larger structures that can potentially receive 
synaptic contacts from many more neurons. Moreover, with 
strong self- avoidance, the surface (space) is almost completely 
covered. This better coverage and larger structure directly results 
from self- avoidance among branches, which therefore grow and 
branch in parallel. The inertial bias alone also covers the surface 
(space) well, but less densely as the total length is considerably less 
when compared to the soma-tropic and self- avoidance examples. 
By means of this analysis, we conclude that sparse coverage of a 
space can be achieved by inertial forces, while dense coverage can 
be created by strong self- avoidance and strong soma-tropism. 

3.2. DENDRITIC MORPHOLOGIES ARE SHAPED BY HOMOTYPIC 
FORCES 

We established a relation between our modeled homotypic forces 
and the resulting morphologies of dendritic arbors. Do these 
traits also exist in real neuronal morphologies? We invert the 
question by generating dendritic morphologies exclusively based 
on homotypic growth biases. If the resulting morphologies 
resemble true morphologies, we can conclude that homotypic 
forces are sufficient to shape dendrites (though they may not do 
so exclusively). We generated three types of virtual neurons based 
on (and subsequently compared to) alpha motor neurons (cat), 
dentate gyrus granule cells (mouse), and Purkinje cells (mouse). 
We took experimentally reconstructed morphologies from the 
NeuroMorpho.org repository (Ascoli, 2006); the motor neurons 
were from Alvarez et al. (1998), the granule cells from Vuksic et al. 
(2008), and Purkinje cells from Martone et al. (2003). We chose 
motor neurons and granule cells are they are often the subject 
of algorithmic generation (Samsonovich and Ascoli, 2003; Marks 
and Burke, 2007b; Torben-Nielsen et al., 2008) and Purkinje cells 
because they have a peculiar and highly specific morphology. We 
adopted a strategy in which we tuned one parameter configu- 
ration to generate a virtual neuron that matches an exemplar 
reconstructed morphology. Then, due to the stochastic nature of 
the algorithm, we generated an arbitrary number of unique mor- 
phologies from this one configuration file using different random 
seeds. 
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First, we investigated whether the motor neuron morphology 
could be dissected into influences of self- referential forces. To 
achieve a proper fit, we found a very strong soma-tropic growth 
bias that falls off rapidly with distance, together with a moderate 
self- avoidance bias (0.08-0.35) and a slightly stronger inertial bias 
(0.36-0.91) was sufficient. Figure 4A shows three reconstructed 
(top row) and three generated, virtual morphologies (bottom 
row). By visual inspection, the virtual morphologies appear plau- 
sible. To visualize the topology, we also constructed dendrograms 
of one reconstructed and one virtual motor neuron in Figure 4B. 
More quantitatively, we plotted the branch order distribution, 



the distribution of the path lengths (from each terminal tip to 
the soma), and the Sholl- intersections for one reconstructed and 
one virtual neuron (Figure 4C). While the distribution of the 
path lengths and the Sholl- intersections are similar, a clear dif- 
ference can be observed between the branch order distribution 
of the virtual and reconstructed neurons. The mismatch in the 
order is a direct consequence of the simplified branching mech- 
anism: we use a slightly modified Galton-Watson process that 
does not mimic the details of the true branch point distribution. 
Similarity between the Sholl- intersections and the distribution of 
path lengths indicates an overall similar pattern of outgrowth: the 





aaa *aa coa aoa uacn izcn usaa 



FIGURE 4 | Alpha motor neurons and their virtual counterparts. 

(A) Reconstructed alpha motor neurons (from Alvarez et al., 1998) on the top 
row and generated virtual neurons on the bottom row. (B) Dendrogram of 
one reconstructed and one generated motor neuron on the left and right, 



respectively. (C) Comparison between reconstructed and generated motor 
neurons based on their measured distribution of the branching order, path 
lengths to terminal tips, and Sholl-intersections. Scale bars in (A) indicate 
250 |xm. 
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same space is occupied and the density of the dendritic trees are 
similar. The fractal dimension [on a branch to branch basis, as 
measured by the Hausdorff fractal dimension (Marks and Burke, 
2007a)] is 1.09 ±0.02 and 1.08 ±0.02 for the reconstructed 
and generated motor neuron in Figure 4A (first column). The 
related contraction (i.e., Euclidean distance over path length) is 
0.73 |xm ± 0.08 and 0.77 |xm ± 0.07 for reconstructed and gen- 
erated morphologies, respectively. Together, fractal dimension 
and contraction are a good proxy for the space-filling proper- 
ties of neurons. Therefore, we can observe that the generated 
morphologies fill the space to the same degree as the real motor 
neurons. Additionally, the bounding space is also set to match 
the space occupied by the exemplar reconstructed neurons and 
therefore the overall dimension of the generated morphologies 
matches well. The path length to the terminal tips is similar 
between experimentally reconstructed (792 |xm ± 196) and gen- 
erated morphologies (836 |xm ± 214). Thus the geometry of the 
generated motor neurons matches the geometry of the exemplar 
reconstructed cells. Because of the good approximation of the 
overall geometry we can conclude that homotypic forces alone 
are sufficient to account for motor neuron geometry. 

Secondly, we set out to generate virtual granule cells. We tuned 
one parameter configuration that matched with GranuleCell-Nr3 
from Vuksic et al. (2008). By generating morphologies from the 



same parameter set but with a different random seed, we could 
generate a set of unique morphologies matching well with the 
population of 20 exemplar cells. In comparison to alpha motor 
neurons, granule cells required a strong soma-tropic growth 
bias with a more gradual fall- off, a considerably stronger self- 
avoidance bias (4.5), and a comparable inertial bias (0.45). Strong 
self avoidance has also been demonstrated experimentally among 
Purkinje cell branches (Fujishima et al., 2012). In addition, the 
flatness parameter was reduced from 1.0 (fully 3-D) to 0.03 
(strongly restricted to 2-D), since we observed that granule cells 
morphologies were partially flattened. A summary of the results 
is illustrated in Figure 5. The visual resemblance between the 
exemplar morphology and the generated morphologies is high. 
The number of terminal tips (14 ±4) and bifurcation points 
(13 ± 4) is similar to the example cell (terminal tips 13 ± 3 and 
bifurcation points 12 ± 3). The total length of the generated 
granule cells, 1024 |xm ± 373 (vs. 2068 |xm ± 241), is slightly less 
for the generated morphology, while path length (205 |xm ± 32) 
(vs. 225 |xm ± 22) is comparable. The fractal dimension is much 
lower than for the motor neurons and amounts to 1 ± 0.02 
for the generated cells and 1.01 ± 0.01 for the exemplar cells. 
From Figure 5 we can see that the distribution of the branching 
order and the Sholl-intersections are similar between the exem- 
plar and generated cells. The path length, on average, is also 
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FIGURE 5 | Dentate gyrus granule cells and their virtual counterparts. 

(A) Exemplar morphology from Vuksic et al. (2008) with its dendrogram. 

(B) Three morphologies generated from one parameter configuration each 



with its associated dendrogram. (C) Comparison between properties of one 
exemplar granule cell (from A) and one generated morphology (first from B). 
Scale bar in (A) and (B) indicate 100 |xm. 
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similar but shows far less variance in the generated cells compared 
to the exemplar. However, based on the overall similarities, the 
shared dimensions, and the similarity in the space-filling prop- 
erty, we can conclude that granule cell morphology is sufficiently 
explained in terms of homotypic forces. 

Lastly, we generated a set of Purkinje cells based on one exem- 
plar cell, namely e4cb2a2 from Martone et al. (2003) (Figure 6). 
Purkinje cell morphology is characterized by a large trunk from 
which smaller branches sprout. Also, the dendritic arbor densely 
covers the plane it occupies, thus requiring that smaller branches 
often grow toward the soma. We adjusted the parameter settings 
to match a generated morphology to that of e4cb2a2, then gener- 
ated a set of different cells from these parameters using different 
random seeds. Specifically, we included a parameter that rendered 
branch diameters following a bifurcation unequal (while main- 
taining their specified Rail's ratio), and the algorithm maintained 
this asymmetry for much of the arbor. To achieve appropriate 
branch orders and good coverage, we also set the front exten- 
sion parameter to a negative value, causing segment lengths to 
vary inversely with radius and thus branching frequency and tor- 
tuosity to slightly increase as the arbor expanded. In comparison 
to alpha motor neurons and granule cells, Purkinje cells required 
a 10-100 times greater branching probability, a moderate soma- 
tropic growth bias with a gradual fall- off, and both the strongest 
self- avoidance bias (5.4) and weakest inertial bias (0.32) among 
the three types. In addition, the flatness parameter was 0.0 



(fully 2-D), since we observed that Purkinje cell morphologies 
were either fully within a plane, or occupied two distinct planes 
(a curious phenomenon our algorithm could not accommodate). 

A quick visual inspection shows great correspondence between 
the exemplar morphology (Figure 6A) and the generated neu- 
rons (Figure 6B). The dendrograms are complex and difficult 
to compare, but the primary feature of a main trunk and a 
higher number of bifurcations combined with a high centrifu- 
gal order are observed. We then compared the complete set 
of four exemplar neurons available from the same lab through 
NeuroMorpho.org. The overall topology of the generated neurons 
compare well to the generated ones [373 =b 84 vs. 354 zb 11 termi- 
nal tips and partition asymmetry, Equation (2), of 0.44 ± 0.02 
vs. 0.53 =b 0.01]. The geometry is also comparable in total length 
(generated 4377 ±918 vs. exemplar 5227 =b 544) and in path 
length (generated 108 =b 10 vs. 126 =b 14), albeit, the generated 
neurons being slightly smaller. This difference in size is in part due 
to the fact that the selected exemplar cell was the smallest of the 
set with a total length of 4461 /xm. We find that the spatial embed- 
ding of the generated neurons as measured by the fractal dimen- 
sion (generated 1.14 zb 0.05 vs. 1.09 =b 0.03) and the contraction 
(generated 0.74 zb 0.04 vs. 0.68 zb 0.02) also matches well with 
the exemplar cells. Curiously, the fractal dimension of the exem- 
plar cells is lower than the generated cells indicating straighter 
branches in the exemplar cells while the contraction indicates 
exactly the opposite (stronger contraction for the exemplar cells). 
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FIGURE 6 | Cerebellar Purkinje cells and their virtual counterparts. 

(A) Exemplar morphology from Martone et al. (2003) with its dendrogram. 

(B) Three morphologies generated from one parameter configuration each 



with its associated dendrogram. (C) Comparison between properties of one 
exemplar Purkinje cell (from A) and one generated morphology (first from B). 
Scale bars in (A) and (B) indicate 100 |xm. 
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FIGURE 7 | Parameters differing across cell types. (A) Force parameters 
plotted on logarithmic scale /-axis for the three cell types: alpha motor 
neurons (OK granule cells (□), and Purkinje cells (A). All parameters plotted 
and superimposed on box plots, with median indicated by a horizontal line, 
the box's edges indicating the 25th and 75th percentiles, the whiskers 
extending to all points not considered outliers, and outliers (>1.5 x the box 
height beyond the box edge) marked by a dot. (B) Spatial gradient 
parameters for the three cell types plotted on a linear scale /-axis. 
(C) Additional relevant parameters for the three cell types plotted on a 
linear scale /-axis. 



A summary of these data is also given in Figure 6C, where the 
match can also be visually evaluated. 

4. DISCUSSION 

In this work we investigated the influence of homotypic forces on 
the morphology of dendrites. We followed a generative synthetic 
approach, in which artificial dendritic morphologies were solely 
determined by homotypic influences. By approximating three 
homotypic forces (inertial, soma-tropic, and self-avoidance) in a 
phenomenological growth model we could generate morpholo- 
gies that have strong resemblance with real dendritic morpholo- 
gies. We conclude that these homotypic forces are sufficient to 
explain dendritic morphology, though they may not do so exclu- 
sively. Consequently, we also argue that an estimate of these 
three homotypic forces should be used in future descriptions of 
dendritic morphology. 

Despite a considerable amount of experimental evidence for 
external cues influencing dendritic morphology, there is little the- 
oretical work on this topic. Moreover, in algorithmic approaches 
to generate whole morphologies, external influences are gener- 
ally not considered. Recent noteworthy exceptions come from 
Samsonovich and Ascoli (2005), Luczak (2006), and Marks and 
Burke (2007a,b) to which we compare our work. Samsonovich 
and Ascoli set to dissect the main direction of growth in hip- 
pocampal pyramidal cell dendrites. They found that the main 
component was a soma-oriented tropism, and not a general 
gradient shared with other cells. Marks and Burke came to a 
similar conclusion and reported that, in motor neurons, the 
direction of the dendrites is both determined by the direc- 
tion of the parent (i.e., what we call "inertial force") and a 
soma- repulsive cue. Moreover, Marks and Burke also hinted at 
the influence of self- avoidance during the development of den- 
drites. An important distinction between first, the algorithms of 
Samsonvich and Ascoli and second, that of Marks and Burke, 
is that each starts with a dendrogram and only then adds a 
geometry. We generate both a topology (i.e., a dendrogram) 
and the spatial embedding (i.e., a geometry) at once, such that 
certain branching events resulted in one branch failing due to 
termination, while the other continues. Luczak used an algo- 
rithm based on diffusion limited aggregation that implicitly 
mimics precisely certain external cues (external gradients and 
self- avoidance). The results of Luczak indicate that a bounding 
box and competition over resources can account for dendritic- 
specific features that resemble self-referential and environmental 
interactions. 

What is the relevance of the presented work? We can assess the 
relevance of the two pillars of our work, the branching rule based 
on a Galton- Watson process and homotypic forces defining the 
geometry. First, a Galton-Watson process is the simplest process 
to generate binary trees and relies on a constant branching proba- 
bility independent of any other parameters in the algorithm. The 
simplicity of this process gives rise to discrepancies between true 
morphologies and the morphologies generated by our algorithm, 
namely the maximum branch order and distribution of bifurca- 
tions. With a pure Galton-Watson process, we would employ a 
constant branching probability (and hence a constant segment 
length independent of other variables and uniform across the 
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whole tree). Several studies have indicated that the branching 
probability is not uniform, and can be dependent on the dis- 
tance from a previous branch point (Nowakowski et al., 1992) or 
from the soma (Burke et al., 1992). Departing from our original 
intent, our algorithm may be modified in the future to generate 
highly realistic morphologies by substituting the Galton- Watson 
process with a more data-driven branching rule, or even a branch- 
ing growth model as developed in van Pelt et al. (2001). The 
modified Galton-Watson process we used here did provide some 
ability to control a non-uniform branching probability, which we 
exploited to approximate more closely the properties of our real 
target neurons, such as their distributions of path lengths. 

Second, homotypic forces were modeled to guide the geomet- 
rical extent of the artificial dendrite morphologies, and varied 
greatly in comparison to other parameters across the different 
cell types we generated (Figure 7). Despite being a phenomeno- 
logical model, the homotypic forces are biologically inspired. 
Stiffness resulting from an inertial force is documented (Condron 
and Zinn, 1997) and biases forward motion. Soma-tropic force 
has been exclusively described in a theoretical setting based on 
a transient gradient of pH induced by a cell's spiking activity 
(Samsonovich and Ascoli, 2005). Self- avoidance is also well doc- 
umented based on observations of homotypic repulsion after 
actual contact (Jan and Jan, 2010; Fujishima et al., 2012) as well 
as from a distance (Lefebvre et al, 2012). In our model, repul- 
sion follows a gradient such that actual contact between branches 
is not required. In addition, Shimono et al. (2010) and Sugimura 
et al. (2007) created detailed mechanistic models of growth at the 
molecular level, which also confirmed that self-repulsive effects 
are important in generating dendritic morphology. In contrast, 
the present study showed that simple models of self- referential 
forces can create realistic and diverse neuronal morphologies 
and did not attempt to capture the proximal physiological or 
molecular causes of these growth processes in a generative algo- 
rithm. Thus, the phenomenological incorporation of inertial, 
soma-tropic, and self- avoidance forces into our morphogenetic 
algorithm are grounded in the literature. 

Of course there are also extrinsic actors (as opposed to self- 
referential or intrinsic forces) involved in shaping dendrites (Scott 
and Luo, 2001; Jan and Jan, 2010). These external influences 
can have a wide variety of sources. Neurons develop in a sub- 
strate populated by both other neurons and other structures such 
as blood vessels, capillaries and glia cells. Especially for opti- 
mal space coverage and tiling, interactions between cells seem 
required, and avoidance of other structures is evident. Moreover, 
neurons grow in a bounded space where the white matter, the pia, 
and even laminar structures can be seen as external influences on 
neurons. At a phenomenological level, these extrinsic actors can 
be described as forces in a (limited) spatial domain (as in Feng 
et al., 2005). Because we model the forces acting on the direction 
of growth as a sampling bias, we can easily extent the model to 
include other external forces in a similar fashion. 

Our work demonstrates that self- referential forces are suf- 
ficient to explain the main morphological traits of adult den- 
drites. Experimental techniques, including genetic manipulation 
(Lefebvre et al, 2012), provide a means to access mechanistic 
causes of self- referential biases. As these techniques evolve, they 
may provide a means to test the predictions represented by our 



choice of parameters for those generated morphologies presented 
here, or for others derived from our algorithm. Our focus was 
not on generating the most realistic virtual neurons possible from 
a set of experimentally reconstructed morphologies. For exam- 
ple, as already discussed, our modified Galton-Watson branching 
process gives rise to a less realistic, nearly uniform branching fre- 
quency throughout the dendritic tree. In a similar vein, other 
specific features of dendritic morphology (e.g., segment length 
and branch angles) are not necessarily uniform across the tree 
(Burke et al., 1992; Nowakowski et al., 1992; Torben-Nielsen et al., 
2008; Langhammer et al., 2010), and instead vary as a function 
of some topological measures of the tree (such as path length, 
branch order, etc.). Therefore, a potential extensions to this work, 
aimed at the generation of realistic morphologies, could include 
other parameters (such as those describing homotypic forces and 
their influence) that vary over the dendrite in this way. 

Recently it was also reported that external guidance cues are 
not required for synaptic patterning and contact formation (Hill 
et al., 2012). While the finding is surprising, it does not speak 
to the need for intrinsic or extrinsic cues during development 
of the morphology itself. Moreover, it is reasonable to think that 
evolution caused a particular trait in a population, giving rise to 
a statistically relevant number of appositions between dendrites 
and axons, to become the basis for regular cell-to-cell wiring and 
communication in neural tissues within a species. This result also 
implies that neurons can, in principle, be generated indepen- 
dently from each other, as long as they share population statistics 
with a measured population of real neurons. In our study we 
managed to dissect morphological traits of dendrites stemming 
from the interaction and influence of three self- referential forces. 
In combination with adequate branching rules, our algorithm 
provides a parsimonious description to generate vast numbers of 
virtual neurons. Moreover, we speculate that morphologies from 
many distinct classes exhibit morphological traits that result from 
self- referential interactions biasing an otherwise random process. 
Tuning the small set of parameters underlying our algorithm is 
more straightforward than tuning an exhaustive set of detailed, 
local parameters (such as the distribution of branching angles, 
which in turn might be conditioned on various other dendritic 
parameters such a diameter and distance from the soma). Thus, 
our insights into the growth of single neurons could be used to 
generate large populations of neurons and synapses for large-scale 
computer simulations of neural tissue. 
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